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Abstract. During the violent relaxation of a self-gravitating system a significant fraction of its mass may be 
ejected. If the time varying gravitational field also breaks spherical symmetry this mass can potentially carry 
angular momentum. Thus starting initial configurations with zero angular momentum can in principle lead to a 
bound virialized system with non-zero angular momentum. We explore here, using numerical simulations, how 
much angular momentum can be generated in a virialized structure in this way, starting from configurations of 
cold particles which are very close to spherically symmetric. For initial configurations in which spherical symmetry 
is broken only by the Poissonian fluctuations associated with the finite particle number N, with N in range 10 3 
to 10 5 , we find that the relaxed structures have standard “spin” parameters A ~ 10 3 , and decreasing slowly 
with N. For slightly ellipsoidal initial conditions, in which the finite-./V fluctuations break the residual reflection 
symmetries, we observe values A ~ 10 -2 , of the same order of magnitude as those reported for elliptical galaxies. 
The net angular momentum vector is typically aligned close to normal to the major semi-axis of the triaxial 
relaxed structure, and also with that of the ejected mass. This simple mechanism may provide an alternative, 
or complement, to “tidal torque theory” for understanding the origin of angular momentum in astrophysical 
structures. 
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1. Introduction 

Observations have shown already several decades ago that 
galaxies of all types generically have significant angular 
momentum. Its origin remains a fascin ating open theoret- 
ical problem (for a review, see e.g. iRomanowskv fe Falll 
(12012 11. Globular clusters have al so been more recentl y 


observed to have net rotation (I Henault- Brunet et al 


20121 : Bellazzini et al. . 2012 : Kacharov et all 2014h . The 
most popular theory to account for angular momentum 
in virialized structures is so-called “tidal torque theory” 
in which the virialized structure gains angular momen¬ 
tum by the action of the torque due to the tid a l field s 
generated by surrounding structures ( Pceblcsl . Il960ll . 
We explore here a distinct mechanism which, despite 
its simplicity, appears not to have been considered 
in the literature, apart from in on e recent study of 
cold spherical collapse by one of us (IWorrakitpoonpon . 
20 lfifl : generation of angular momentum by ejection of 


matter in violent relaxation. Indeed violent relaxation 
of self-gravitating structures is characterized by very 
large amplitude variations of the mean field potential of 
a structure, which can give enough energy to particles 
to escape even if all the mass is initially bound. The 
amount of mass ejected depends strongly on the initial 
conditions, varying from zero to as much as 40% for 
highly unif orm and com pletel y col d initial conditions 
( Joyce. Marcos fc Svlos Labinil . 2009lf1 . If, in addition, 
the mass distribution is not spherically symmetric during 
the violent phase of the relaxation, this ejected mass 
would be expected to carry at least some angular mo¬ 
mentum — and “generate” in the remaining bound mass 
an (equal and opposite) amount of angular momentum. 
In other words the two components of the mass — the 
part that is ejected, and the other part which remains 
bound — can exert a net torque on one another during 
the violent relaxation leading to a net angular momentum 


1 Mergers of virialized structures can also lead to ejection of 
a significant mass, i n particular when th e two structures have 
very different mass (ICarucci et all [2014 1. 
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for both of them. Further, given that violent relaxation 
starting from a wide range of cold initial conditions is 
often characterized by a very strong breaking of spherical 
symmetry even when the initial conditions are not 
leadin g in parti cular to t riaxial relaxed structures (e.g. 
Merritt & Aguilar (tl985l l: Barnes. Lanzel fc William; 


( 2009h : _ Benhaiem &; Svlos Labinil ( 20151 , 

Svlos Labini. Benhaiem fe Joyce ( 20151) 1 - - one might 


expect the effect to be far from negligible. 

We explore here, using numerical simulations, how 
much angular momentum can be generated in a virial- 
ized structure by this mechanism. More specifically we 
consider an isolated initially cold distribution of matter in 
open boundary conditions, and without expansion, and a 
range of initial spatial distributions which are very close 
to spherically symmetric: (i) particles distributed around 
a center following a mean density profile which decays as 
a power law of the radial distance, or (ii) particles dis¬ 
tributed with uniform mean density inside an ellipsoidal 
region. These are initial conditions which have been ex¬ 
tensively studied in the literature (see references below) 
to investigate the processes involved in the formation of 
galaxies and other astrophysical structures. We note that, 
in a cosmological context, the simulations can be taken to 
represent the evolution, in physical coordinates, of a single 
isolated overdensity with the chosen initial profile. On the 
other hand the relation of such simulations to the more 
general case of the evolution an overdensity in an expand¬ 
ing universe (which cannot necessarily be well approxi¬ 
mated as isolated) is a more complex issue, and will be 
discussed in detail in a separate forthcoming publication. 
It involves several additional parameters which can play 
a role in the dynamics (e.g. how precisely the background 
density is modelled) and additional numerical issues (no¬ 
tably concerning the control of numerical accuracy wit h¬ 
o ut energy conserva tion, see Joyce & Svlos Labinil (12012 1: 
Svlos Labini ( 2012h J. 

In the present context the fluctuations breaking spher¬ 
ical symmetry in the initial conditions are evidently of 
central importance for the phenomenon we are studying. 
In the first case spherical symmetry is broken only by the 
finite particle number fluctuations, while in the second 
case these fluctuations break the residual reflection sym¬ 
metries.We find that almost all these initial conditions 
except where there is negligible mass ejection indeed 
lead to a measurable net angular momentum in the re¬ 
laxed virialized structure. The magnitude of this angular 
momentum is larger by about an order of magnitude in 
the latter class, and in many cases is sufficiently large 
to suggest that the mechanism could potentially account 
for angular momenta observed in astrophysical structures. 
Indeed, the typical meas ured values of th e spin parame- 
ters of galaxies (see e.g. iHernandez et al.l ( 20071') ') are of 
the same order of magnitude as those we find in our sim¬ 
ulations. 

The paper is organized as follows. We first present the 
details of the numerical simulations and initial conditions 
in Sect|2j In the following section we then present our re¬ 


sults, first describing the relevant features of the evolution 
qualitatively and then giving the quantitative results. In 
Sect|T]we summarize our results and conclude. 


2. Numerical simulations 

2.1. Initial conditions 

In detail the initial conditions of which we study the evo¬ 
lution are the following: 


N particles distributed randomly inside a sphere of ra¬ 
dius Rq, following the radius-dependent density pro¬ 
file p(r) (x r~ a where r is the radial distance from 
the cent re and a is a constant. We will refer to 
these as “spherical initial conditions”. The particles 
have vanishing initial velocities. We report here re¬ 
sults for the range of a with 0 < a < 2. We re¬ 
strict to this range as m a previous study dSvlos Labini 


2012 : Svlos Labini. Benhaiem fe Joyce . l2015f) we have 
observed that, for a > 2 , there is negligible (-C 1 %) 
mass ejection. As shown in this same study, the evo¬ 
lution leads, except for a very close to zero, to virial 
equilibria which are very non-spherically symmetric, 
and typically triaxial. We have varied the number of 
particles N from 10 2 3 * to 10 5 . 

N particles are distributed randomly in a prolate el¬ 
lipsoidal region. We define <23 to be the largest semi¬ 
principal axis and ai = 02 the smaller ones. The three 
eigenvalues of the inertia tensor are simply 


A., 


= 


4 ) 


(i) 


where M is the total mass of the system, and i ^ j ^ k 
and i,j,k = 1,.., 3: from the definition of the semi¬ 
principal axes we have Ai > A 2 > A 3 . Note that the 
principal axis corresponding to the eigenvalue Ai is 
oriented in the direction of the shortest semi-principal 
axis ai, while the principal axis associated with A 3 
is in the direction of the longest one 03 . The particle 
number N spans again the range from 10 3 to 10 5 . 

The shape at any time may then be characterized by 
the parameter d 


i{t) 


Ai (t) 
A 3 (t) 


( 2 ) 


We will refer to these as “ellipsoidal initial conditions”, 
and report here results for value of 1 in the range from 
t(0) = 0.01 to t(0) = 0.25. The mass ejection and am¬ 
plification of the spherical symmetry breaking during 
the evolut ion from t hese i nitial condition s has been 
studied in Benhaiem fc Svlos Labinil (2015:). 


2 Note that these parameters are generally defined as a func¬ 

tion of the semi-principal axes ffli, 02 , 03 rather than as a func¬ 

tion of the eigenvalues (see EqQ. However for small deforma¬ 

tions of a perfect sphere, which is the case we consider here, 
these definitions are almost equivalent. The advantage of this 
definition of the parameters using the eigenvalues is that it can 
be used to characterize any distribution. 
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2.2. Numerical simulations 


We have used the N-body code Gadget ( Springe! . 12005 ). 
All results presented here are for simulations in which en¬ 
ergy is conserved to within one tenth of a percent over the 
time scale evolved, with maxi mal deviations at any time 
of less than half a percent (see Benhaie m fe Svlos Labini 
( 20151) : ISvlos Labini. Benhaiem fc .Tovce ( 2015 ) for more 
details). This accuracy has been attained using values of 
the essential numerical parameters in the Gadget code 
[0.025 for the rj parameter controlling the time-step, and 
a force accuracy of of = 0.001] in the range of typi¬ 
cally used values for this code. In the specific cases of 
spherical initial conditions with a = 0, which is sin¬ 
gular in the limit N —> oo, the treatment is as de- 
tailed in Joyce. Marcos fc Svlos Labini ( 20091) (see also 
Iworrakitpoonpon ( 2015 )). We have also performed exten¬ 
sive tests of the effect of varying the force smoothing pa¬ 
rameter e, and found that we obtain very stable results 
provided it is significantly smaller than the minimal size 
reached by the whole structure during collapse H. We will 
discuss the dependence on particle number N of our re¬ 
sults below. 

As noted above, our simulations are performed with 
open boundary conditions and in a non-expanding back¬ 
ground, but can be taken to represent well the evolution of 
an isolated overdensity in an expanding universe, which, at 
the initial time, is of high density compared to the mean 
mass density of the universe and at rest, in physical or 
comoving coordinates. In this respect we could perform 
the simulations here using the expanding universe version 
of the code and periodic boundary conditions. The size 
of the periodic box relative to that of the initial sphere 
would then fix the relative overdensity represented by the 
sphere at the initial time (and thus also the time scale for 
virialization compared to the Hubble time). Such a simu- 
lation gives, as discussed in detail and stud i ed at l ength in 
Joyce fc Svlos Labin i ( 2012 ): Svlos Labinil ( 20121) . identi¬ 


cal results in physical coordinates, modulo finite size cor¬ 
rections suppressed by the ratio of the size of the struc¬ 
ture to the size of the periodic box, and possible numer¬ 
ical effects (incl uding the effect o f s moot hing). I ndeed, 
as dis cussed in i Joyce fc Svlos Labini ( 2012 ): Svlos Labini 
( 20121 ) . the differences between such simulations and those 
in open non-expanding space provides information about 
such numerical effects, and indeed a tool to control bet¬ 
ter expanding universe simulations. Our “direct” simula¬ 
tions in physical coordinates are preferable because they 
are simpler and, notably, far easier to control for accuracy 
(via energy conservation). 

Given that we are interested here in the angular mo¬ 
mentum, we also test our simulations for its overall con¬ 
servation and, as detailed below, compare the measured 
angular momentum of the final bound structure with the 


3 More specifically we have found very stable results for e in 
the range R™ 1 ™ /e e [10, 200] where R™ tn is the minimal value 
of the gravitational radius (defined further below) reached dur¬ 
ing the collapse. 


numerical error. We will also check systematically for the 
accuracy of conservation of the total linear momentum. 


3. Results 

3.1. Qualitative description of mechanism 

Numerical studies of evolution from some of these initial 
conditions, or very similar ones with small but non-zero 
virial ratios, have been re porte d extensively in th e 


literature (see 


Tigj 


Henon 


Aarseth. Lin fc Panaloizou 


( 19731) 

A 


van Albada (1982D: 


Theis & Spurzcml 


(1999 _ Boilv. Athanassoula fc Krounal (20021): 


unt oo ivioumu| _ 

Joyce Marcos fc Svlos Labinil (l2009h : Svlos Labinil 


(12012 . _ 20131): iBenhaiem fc Svlos Labinil ( 20151) : 

Svlos Labini. Benhaiem fc Jovcel ( 2015 )), with an 
emphasis in particular on the study of the shape and 
profile of the virialized structure. Because the system is 
initially cold it undergoes in all cases a strong collapse, 
in a time of order 1 / y/Gpo where po is the initial mean 
density, followed by a re-expansion which leads rapidly to 
virialization of most of the mass. The “degree of violence” 
of the collapse — which can be characterized roughly by 
the maximal contraction the system undergoes — varies 
with the initial condition. The most extreme case is the 
spherical case (a = 0) in which the collapse is singular in 
the limit N —> 00 . 

While both early studies, and most studies since then, 
have focused on the properties of the virialized struc¬ 
ture (e.g. its density profile, shape, and velocity dis¬ 
tributions) the phenomenon of mass (and energy) ejec¬ 
tion has been considered in detail only more recently: 
for the case of cold unif orm spherical conditions (i.e. 
the ca se a = 0 here) in Joyce. Marcos fc Svlos Labinil 
(120091) . for the case a = 0 with non-zero initial veloci¬ 
ties in Svlos La binil (l20 12h and for el lipsoi dal initial con¬ 
ditions in Benhaiem fc Svlos Labinil ( 2015 ) Q. Essentially 
the amount of mass ejected depends on the degree of vi¬ 
olence of the relaxation, which can be quantified roughly 
by the maximal contraction attained by the system dur¬ 
ing the collapse phase. The reason for this correlation 
between the ejection of mass and the violence of col¬ 
lapse becomes easy t o understand when the mechan i sm for 
ejection is stu died (Joyce, Marc os fc Svlo s Labinil . 2009; 


Carucci et all 20l4 Benhaie m fc Svlos Labini 


2015 


Svlos Labini. Benhaiem fc Jovcel 2015 ): The ejected par¬ 


ticles are essentially those of which the fall times to the 
center of the structure during collapse are longest. As a 
consequence they pass through the center of the whole 
structure when it has begun re-expanding, which means 
that they travel into a deeper potential than the one that 
they travel out of. As a result they pick up an energy 
“kick”. The magnitude of this energy gain is directly re¬ 
lated to the strength of the potential at this time, which 
depends essentially on the extent of the contraction. 


Mass ejecti on in mergers o f two v i rialized st r uctur es is stud¬ 
ied in detail in lCarucci et al.l (120141 ): ISamsinel (120151 ). 
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For the generation of angular momentum, the sec¬ 
ond crucial ingredient, in addition to mass ejection, is 
the breaking of spherical symmetry. Indeed mass ejec¬ 
tion can occur in a purely radial gravitational field, 
but will not then generate angular momentum. That 
cold initial conditions in the class we are considering 
lead to virialized states that are ver y far from spheri¬ 
cally s ymmetric was observed first by iMerritt fe Aguilan 
( 1985ll. and extensive l y stud i ed in the literature (se e 
e.g. A guilarj^^Mgrritt] ( 1990h : Theis fc Sourzem (19991: 


Boilv fc Athanassoulal (teOOfill : iBarnes. Lanzel fc Williams! 

(l2009h j.'This instability of initially spherical systems to 
relaxation to non-spherical viral equilibria is usually re¬ 
ferred to as “radial orbit instability”, because of its ap¬ 
parent link to an instability of equilibri um system s with 
purely radial orbits o riginally proved bv lAntonovI (|19611) : 


Fridman et al 


. In two recent papers we have stud¬ 
ied in detail the development of this asymmetry during 
the collapse phase for both the classes of initial con¬ 
ditions we study in this paper (s p herica l initial condi¬ 
tions in Svlos Labini. Benhaiem fe Joyce rt2015lh. and el ¬ 
lipsoid al initial conditions in Benhaiem fe Svlos Labinll 
( 2015IB . These studies reveal that the same process in¬ 
volved in the energy ejection leading to mass ejection plays 
a crucial role, for many initial conditions, in amplifying the 
symmetry breaking. Indeed, amongst the late arriving par¬ 
ticles, there is also a spread in arrival times as a function of 
angle, which leads to an energy injection, and thus a spa¬ 
tial distribution of mass, which is also a function of angle. 
The effect is maximal in enhancing the final asymmetry in 
the range a £ [0.5,1.5]. It is, on the other hand, relatively 
suppressed in the limit a = 0 because the particles’ fall 
times, and consequently their energy changes, are much 
less strongly correlated with their initial radial position in 
this case, leading to a “washing out” of the effect. The de¬ 
velopment of the asymmetry in the mass distribution, and 
thus in the gravitational potential, is the combined result 
of the growth of the initial finite N fluctuations break¬ 
ing spherical symmetry thr ough gravitational i nstab ility 
- known in this case as the Lin. Mestel fc Shu ( 1965ll in¬ 
stability — which is then amplified by the energy ejection 
to particles which also leads to mass ejection. 

A schematic representation of the mechanism is given 
in Figure [T] Any given realization of our initial conditions 
has a longest semi-principal axis, along which particles are 
on average further from the orthogonal plane. In the case 
of the spherical initial conditions the non-zero effective el- 
lipticity is a finite N effect, while in the ellipsoidal initial 
conditions it is dialed by the parameter (.(0). In the first 
phase of collapse this initial a symmetry is amplif i ed by 
the (gravitational) instability of Lin, Mestel fc Shul ( 19651) 
with the collapse occurring latest along the longest axis. 
Particles arriving from the corresponding directions pass 
through the center as the structure is already re-expanding 
in the other directions, leading to a greater energy injec¬ 
tion to them. As the potential they are traveling is not 
spherically symmetric the particles also gain traverse ve¬ 
locities and thus angular momenta with respect to the 




Fig. 1 . Schematic representation of the generation of an¬ 
gular momentum in the ejected particles. The upper panel 
shows schematically symmetry breaking in the initial 
cloud. When particles are ejected both their radial and 
traverse components will generally have some dependence 
on angle, leading to a net angular momentum being car¬ 
ried away by the ejected particles. 


center. As there are also fluctuations breaking rotational 
symmetry, which grow in the course of the collapse, both 
the radial and transverse velocities will vary as a function 
of direction, and, for example, those ejected in opposite 
directions will have slightly different angular momentum. 
The ejected outgoing particles can then carry net non-zero 
angular momentum leaving behind the opposite angular 
momentum in the particles which virialize in the bound 
structure. 

3.2. Measurement of angular momentum 

We now focus on the evolution of the angular momentum 
during the relaxation. More specifically we decompose the 
total angular momentum Lt at any time into two compo¬ 
nents: 

L T = L b + L f (3) 

where L b (Lf) is the total angular momentum of the 
“bound” (“free”) particles, i.e., which have negative (pos¬ 
itive) energy at the given time. In practice particles’ en¬ 
ergies almost never change sign more than once, i.e., the 
asymptotically ejected particles are those which are tagged 
as “free” when their energy becomes positive. Further we 
consider the decomposition of L b as 

L b = Lf m + LI (4) 

where Lf m = M b R b x F(, is the angular momentum of 
the center of mass, located at R b and moving with ve¬ 
locity V b , of the bound particles, and L% is the angular 
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momentum of the bound particles with respect to their 
center of mass. It is the latter which interests us here pri¬ 
marily, but we also monitor and Lt, measured with 

respect to the fixed origin at the center of mass of the 
initial configuration, to assess the accuracy of the simu¬ 
lation for which the total angular momentum should be 
conserved (and equal to zero). 

We will measure time in units of 

r ‘ = \fi Ik (5) 


where p 0 is the total initial mass density, i.e., the total 
mass divided by the volume of the initial sphere or ellip¬ 
soid. It corresponds to the time for a particle initially at 
the outer periphery of a cold sphere with this mass den¬ 
sity to fall to its center in the continuum approximation, 
(i.e. taking TV —> oo and keeping the initial mass density 
profile fixed) and without shell crossing. 

We will measure angular momentum in the natural 
units given by 


L 0 — 


GM 5 / 2 


( 6 ) 


where M is the total (initial) mass of the system, and 
Eq its total (initial) energy (equal to its potential en¬ 
ergy). Further, wc will report results for the angular mo¬ 
mentum of the bound mass given in terms of the so - 
cal led “spi n paramete r” A as defined bv IPeeblcs ( 196911 : 


Kncb c fc P ower (2008): 


A = 



(7) 


where L b Q is 


gmI ' 2 
~ vW 


( 8 ) 


where Mb is the bound mass and Eb (M4) the total (grav¬ 
itational potential) energy of this mass (with respect to 
its center of masspj. 


3.3. Results for a chosen initial condition 

We present results first for the case of spherical initial 
conditions with a = 1, and N = 10 5 particles. We consider 
this case as its evolution is typical of what we observe 
for all our different initial conditions. The panels of Figj2] 
show for this case, as a function of time, (i) the fraction f + 
of the initial mass in free particles, and the gravitational 
radius defined as R g (t) = (ii) the flattening ratio 

(-80 as defined by ©> with the subscript indicating that the 
inertia tensor is calculated using the particles with energy 

5 Another c ommo nl y used norma lization for the spin is that 
introduced in iBullock et al.l (l200lh . defined by A' = \L^\/L\ 

with L\ = sj2M^GRb with Rb = — where Wt, is the 

potential energy. For a virialized structure, with Eb = W4/2, 
A' = \/5/3X. 



Fig. 2. Evolution from spherical initial condition with a = 
1, and N = 10 5 : (top left) the fraction of particles with 
positive energy / + (full black line) and the gravitational 
radius R g (dotted red line) in units of the initial sphere 
radius Rq ; (top right) the flatness ratio tso for the 80% 
most bounded particles; (bottom left) the total angular 
momentum Lt (full black line) with respect to the initial 
center of mass, the angular momentum of the center of 
mass of the bound particles (red dashed line) and 

the angular momentum of the bound particles (blue 
dashed-dotted line) with respect to their center of mass 
(all in units of L 0 defined in Eq. [Gj ; (bottom right) the 
spin parameter A. 


less than 20 % of the energy of the most bound particle, 
(iii) the magnitude of the angular momenta | Lt\, |X,g orrl |, 
\L^\ (on a logarithmic scale), and (iv) the spin parameter 
A. 


The first plot illustrates clearly, as described above, 
that at r c « 1, the system reaches its maximal contrac¬ 
tion, and the strong energy injection which leads to par¬ 
ticle ejection occurs in a very short interval around this 
time, corresponding approximately to the dispersion in 
fall times of the mass. Likewise we see in the second plot, 
showing the behavior of that the breaking of spher¬ 
ical symmetry grows monotonically from the beginning 
and is then strongly amplified around the maximal col- 
l apse by the energy i njection (a s desc ribed in detail in 
Svlos Labini. Benhaiem fc Jovcel ( 2015t B. 


As anticipated, as can be seen from the lower two plots, 
the bound mass indeed goes toward a stationary state 
with non-zero angular momentum (relative to its center of 
mass). This final angular momentum of the bound mass, 
albeit small in the natural units, is clearly not numeri¬ 
cal in origin: it is about three orders larger than the fi- 
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nal total angular momentum which gives a measure of 
the violation of angular momentum conservation due to 
numerical error. Further we have performed several tests 
varying the essential numerical parameters and found our 
results to be stable. On the other hand the total angular 
momentum L b of the bound mass (with respect to the 
initial center of mass of the whole system) is dominated 
by its translational motion: the bound mass “recoils” with 
a net linear momentum compensating that of the ejected 
mass, along an axis which is off-set from the initial center 
of mass. We will give some further details below on this 
ejected linear momentum, which is measured numerically 
to a precision of order 10 -6 (i.e. the measured momen¬ 
tum change of these particles is fO 6 times larger than the 
change in the total momentum due to numerical error). 


3.4. Generation of angular momentum during collapse 

It is interesting to follow the dynamics leading to the gen¬ 
eration of the ejection of net angular momentum in the col¬ 
lapse phase. During this phase the mass which gets ejected 
has not yet undergone the energy boost which leads to its 
ejection, but the gravitational potential has already devel¬ 
oped very significant asymmetry which leads the particles 
to have non-radial velocities, and thus non-zero angular 
momentum (while the total angular momentum remains 
zero). 

Shown in Fig|3] is, again for the case a = 1 and 
N = I0 5 , the measured dispersion of the transverse ve¬ 
locities in radial shells, as a function of radius, and at 
different times during the evolution. The velocity is nor¬ 
malized in units of vo = \JGM/R 0 . We observe that the 
transverse velocities are already, at t = 0.5r c comparable 
with, and, at t = 0.75r c , even larger than their magni¬ 
tude in the final virialized structure. Note that these are 
the transverse velocities produced by the growth of the 
initial density fluctuations breaking spherically symmetry 
through gravitational instability, well before the energy 
injection to particles arriving around t ss r c which leads 
to the mass (and angular momentum) ejection. At longer 
time scales, i.e. t > t c , the outer shell, corresponding to 
very weakly bound particles with energy very close to zero, 
is still expanding as these particles are still traveling out¬ 
ward on their large orbits. On the other hand the inner 
shell has reached a stationary state, as can be seen by 
comparing the number density profiles at t = 2.5r c and at 
t = 5 T c . 



r/R 


o 


Fig. 3. Average in radial shells of the square of the trans¬ 
verse component of particle velocities (in units in which 
vo = \JGM/Ro is unity), as a function of radius, at each 
of the indicated times, for a simulation with N = 10 5 par¬ 
ticles of spherical initial conditions with a = 1. Note that 
already well before the system reaches it maximal con¬ 
traction, at t. ss 1, the Poissonian fluctuations breaking 
spherical symmetry have been amplified to produce tan¬ 
gential velocities larger than in the final virialized state. 


Shown in the upper panel of Figure [5] is the tempo¬ 
ral evolution of the average value of the modulus of the 
particle angular momentum, for the bound particles (with 
respect to their center of mass, in units of Lq). In line with 
what would be expected from the plots of the transverse 
velocities, we observe a constant growth during the col¬ 
lapse phase, followed by a very significant boost around 
the time of maximal contraction, when the energy ejection 
occurs. Thus we see the angular momentum generated by 
the collapse receives a very significant boost in the final 
phase of the collapse. This is further quantified by the 
lower panel of Figure [5] which shows the evolution, as a 
function of time, of 


l p 

l b 


LI 

Ei W 


(9) 


In FigUis shown the distribution P(|f|) of the absolute 
values of the particle angular momentum £ with respect to 
the center of mass of bound particles, at a few different in¬ 
dicated times (before and after the collapse). We observe, 
in line with the results on the transverse velocities, that 
significant angular momentum is generated already before 
the completion of the collapse phase. Further for t > r c 
P{\l\) rapidly approaches the distribution characterizing 
the final virialized state. 


i.e. the modulus of the angular momentum of bound par¬ 
ticles, normalized to the sum of the moduli of the angu¬ 
lar momenta of the same particles. The dashed horizontal 
line corresponds to the value yjN + /iV_ re f + /\/N , which 
is the amplitude of the bound angular momentum which 
would be expected if the ejection operated as a simple ran¬ 
dom removal of particles without any modification of their 
angular momentum. We observe that the final value is in 
fact at least several times larger. Thus, during the final 
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Fig. 4. Distribution of the absolute values of the particle 
angular momentum t (in units of Lq) with respect to the 
center of mass of bound particles P{ |£|), for a = 1 and 
TV = 10 5 at different times before and after the collapse. 


phase of collapse in which the large energy changes giv¬ 
ing rise to ejection occur, the ejected and bound particles 
exert a torque on one another which very significantly am¬ 
plifies their final (equal and opposite) angular momenta. 

3.5. Results for different initial conditions (fixed TV ) 

In the above we have focused, for simplicity, on the sin¬ 
gle case of spherical initial conditions with a = 1, and we 
have shown results for TV = 10 5 . The qualitative behaviors 
we have discussed in this case, for what concerns the gen¬ 
eration of angular momentum, are shared by all our other 
initial conditions. Quantitatively there is a variation of the 
final angular momentum, which depends notably on the 
details of the mass ejection and the degree of symmetry 
breaking during the collapse. The precision of the simula¬ 
tions also varies from case to case — simulations for the 
spherical initial conditions with a = 0 are particularly del¬ 
icate as the collapse is the most extreme in this case, with 
a singular behavior in the limit TV —> oo. Correspondingly 
in this case the measured angular momentum is only a 
few times larger that the numerical error in the conserva¬ 
tion of the total angular momentum, while it is, as shown 
above (see Fig. ED almost three orders of magnitude for 
the case a = 1. 

Fig [G] shows the “final” values of /+ (top panel), tgo 
(third panel) and A (bottom panel), i.e. at t = 5r c , for 
the spherical initial conditions as a function of a. In addi¬ 
tion, in the second panel, is shown also /g, the “fraction of 
ejected energy”, defined by f E = {E b -E 0 )/E 0 = -E p /E 0 


Fig. 5. For the case a = 1, and TV = 10 5 , the upper 
panel shows the evolution of the average of the modulus 
of the angular momenta of bound particles (with respect 
to their center of mass) as a function of time. The lower 
panel shows the quantity in Eq[9l and the dashed line the 
estimated value expected if the ejection were a random 
sampling of the total mass. 


where E p is (by energy conservation) the total ejected en¬ 
ergy in the frame in which the particles are initially at 
rest, i.e. the sum of the kinetic energy of the center of 
mass of the bound mass and the kinetic energy of the 
ejected particles. FigJT] shows the same four quantities for 
the ellipsoidal initial conditions, as a function of i(0). 

Comparing these two figures, the most striking result 
is that there is an order of magnitude difference in the 
final angular momentum (as characterized by the spin pa¬ 
rameter A) for most of the spherical initial conditions with 
a ^ 0, and the ellipsoidal initial conditions with t(0) 0, 

with a sharp interpolation between the two classes around 
the case a = 0 (which indeed is also the limit i(0) —> 0 
of the ellipsoidal initial conditions). In both figures we 
show also in the lower panel the final angular momentum 
normalized to Lq. Comparing with A, we see that around 
a = 0, where the relaxation is most violent, there is a sig¬ 
nificant difference. This arises simply from the fact that 
there is a considerable ejected mass, and energy, in this 
case: indeed the relative amplitude of the two quantities 
can be expressed as 


A = (K \ 5/2 = ( x + fe ) 1/2 

L\ \M b ) 'vlEoiy (1-/+) 5 / 2 ' 


( 10 ) 


Thus the final spin includes an amplification which comes 
from the change in the characteristic scale of angular mo¬ 
mentum due to the mass and energy ejection (which leads 
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Fig. 6. Final values of /+, iso and A for the spherical initial 
conditions as a function of a, for N = 10 4 . The points are 
the average values for N s = 20 realizations for each a ^ 0, 
and for 47 realizations in the case a = 0; the indicated 
error bar is the corresponding standard deviation a (for A 
we show the error on the mean a/y/N s ). The lower panel 
also shows (red points, dashed) the angular momentum of 
the bound mass normalized to Lq. For the smallest values 
of a the angular momentum generated actually decreases, 
but, because of the increasing mass (and energy) ejection 
the characteristic angular momentum of the final bound 
mass decreases even faster. 


to a more bound, but less massive, structure than the ini¬ 
tial condition). 

It is clear, however, from the two plots that this am¬ 
plification from the normalization is not the main factor 
explaining the much larger final angular momenta from 
the ellipsoidal initial conditions. Rather the main effect 
amplifying the spin in the ellipsoidal initial conditions 
compared to the spherical ones comes clearly from the 
essential difference in the initial conditions: the stronger 
initial breaking of the spherical symmetry. This initial 
asymmetry both leads to a more anisotropic collapse com¬ 
pared to that in the spherical case) and, at the same 
time, the ejection of a comparable amount of mass to 
that in t he s mall a spherica l models. A s described in 
detail in Benhaiem fc Svlos Labinil ( 2015 ). the particles 
on the stretched axis of the initial condition fall slightly 
later, and receive a large energy boost. Note that it is 
only for tso(0) < 0.15 that the final i(t) is linearly pro¬ 
portional to the initial one. Instead, for i-so(0) > 0.15 , 
substructures form during the collapse, leading to a sub¬ 
stantial difference i n the s hape o f the virialized structure 
Benhaiem & Svlos Labinil ( 20ld l. Similarly the fraction of 


Fig. 7. Final values of /+, iso and A for ellipsoidal initial 
conditions, as a function of i(0). The points are averages 
for 20 realizations, and the error bar the corresponding 
standard deviation (for A we show the error on the mean). 
As in the previous figure the lower panel displays also (red 
points, dashed line) the final angular momentum normal¬ 
ized to Lq). 


ejected particles decreases linearly with i 8 O ( 0 ) only for 

i 8 O (0) < 0.15. 

Besides the angular momentum itself, the generation 
of a net relative motion of the center of masses of the final 
bound and ejected mass is also a direct measure of sym¬ 
metry breaking. We would thus expect that we might find 
a similar trend in the final linear momentum of the ejected 
(or bound) particles. Shown in Fig. [5] is a plot, for all our 
simulations with N = 10 4 particles, of both A and /i, the 
modulus of the “final” linear momentum of the ejected 
particles (equal and opposite to that of the bound mass 
up to a numerical precision which is smaller by several 
orders of magnitude in all cases.) We observe that there 
is indeed a clear correlation between the two quantities. 
Further it appears highly consistent with a roughly linear 
relation, as one might expect, given that they are both 
indirect measures of an initially small spherical symmetry 
breaking. 

We have studied finally also the scale dependence of 


the normalized spin parameter defined as ( Bulloc k et al 


2001 ) 


A'(r) - 


L b(r) 


v / 2GM(r) 3 r 


( 11 ) 


where (r) is the angular momentum of bound particles 
with distance < r from their center of mass and M(r) is 
the mass enclosed in a sphere of radius r. Results for five 
different sets of simulations with a = 0,0.5,1,1.5, 2 and 
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Fig. 8. Final spin parameter A, and the modulus of final 
total linear momentum of the bound particles p , in all 
our (322) realizations of initial conditions (both spherical 
and ellipsoidal) with N = 10 4 . p is normalized in units in 
which GM 3 /R 0 = 1. The dashed line is A oc p. 


N = 10 5 are shown in Fig© The profiles are similar in 
all cases, but their amplitude depends on a as highlighted 
in Fig© The fact that A'(r) is larger for a = 0 than for 
a > 0 is related to the fact orbits are slightly more radial 
in the former case. 


3.6. Correlation of angular momentum with final 
spatial distribution 


As ana l yzed i n deta il in iBenhaiem fc Svlos Labinil (120151 ); 
Svlos Labini, Benhaiem fc Joyce ( 20151) there is, for both 
classes of initial conditions, a correlation between the di¬ 
rection in which mass is ejected and the longest axis of 
both the initial condition and the final bound mass. These 
latter two are strongly correlated because it is the growth 
of the initial fluctuations breaking spherical symmetry, 
amplified by the mass ejection, which leads to the final 
spatial asymmetry. The fact that the ejection of mass is 
amplified along these same directions is because they are 
the particles which fall latest and pick up as a result a large 
energy kick. Given that we have seen that most of the final 
angular momentum is generated at the time of ejection, 
we would expect that it would preferably be aligned or¬ 
thogonal to the preferential axis for ejection. 

Shown in Figs. flOl and fill are histograms, for 20 real¬ 
izations of each indicated initial condition, of the modulus 
of the cosine of the angle between the final angular mo¬ 
mentum of the bound mass, L^ 1 and the eigenvector corre¬ 
sponding to the longest axis of the final bound mass (blank 


Fig. 9. Spin parameter A '(r) (see Eg fill), averaged over 
20 realizations, for a = 0,0.5,1,1.5, 2 and N = 10 5 . The 
distance is normalized to the structure size R m ax- 


histograms), and the final ejected mass (filled histograms) 
For the ellipsoidal initial conditions - which very strongly 
single out an initial preferred axis which remains strongly 
correlated with the final longest axis - the anticipated 
correlation is very clearly present. In the spherical initial 
conditions, the correlation is weaker but nevertheless vis¬ 
ible, except perhaps in t he case a — 0. I ndeed in this lat- 
ter ca se, as discussed in ISvlos Labini. Benhaiem fc Joyce 
(1201511 . the symmetry breaking in the final state is in fact 
much weaker, and the correlation between the final and 
initial asymmetry likewise. This is a result of the fact that 
in this case, in which the fall time of all particles is the 
same in the limit N —> oo, there is a much weaker corre¬ 
lation between the initial radial position of a particle and 
the energy change it undergoes in the violent collapse. 


3.7. Dependence on N 

The N dependence of the angular mo mentu m generated 
by ma ss ejection has been explored in IWorrakitpoonoon 
(12015 1 for the case a = 0, and evidence for a monotonic 
decrease ~ 7V -/3 , with a best fit around /3 = 1/3, was 
found. Inspection of the results we have reported above, 
for simulations with N = 10 4 and N = 10 5 , are consis¬ 
tent with this finding for this case, while for initial condi¬ 
tions with a > 0, an apparent decrease of the average L b p 
as N grows is also observed. To explore this further, we 
have performed a more detailed study of the case a = 1, 
performing a larger number of simulations for a greater 
number of values of N in the range 10 3 to 10 5 . We have in 
this case also evolved each simulation for a longer time to 
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leas 0 | leas 0 





Fig. 10. Histogram of the cosine of the angle, well after 
the collapse, between and the eigenvector correspond¬ 
ing to the longest principal axis of the final bound mass 
(blank histograms), and of the ejected mass (filled his¬ 
tograms) , for 20 realizations of spherical initial conditions 
for different a indicated in the caption. 


leas 0 


law 0 




leas 0 | I cos 0 


Fig. 11. Pdf of the cosine of the angle, well after the col¬ 
lapse, between L £ and the eigenvector corresponding to 
the longest principal axis of the final bound (blank his¬ 
tograms), and of the ejected mass (filled histograms), for 
20 realizations of ellipsoidal initial conditions for different 
c( 0 ) indicated in the caption. 


Fig. 12. Logarithm of the final angular momentum of 
bound virialized mass L v h (in units of Lq) as a function 
of particle number N, for spherical initial conditions with 
a = 1. The point for the cases N = 10 5 correspond to 
the results reported above, at t = 5, while for the other 
N there is data at three times, extending to t ~ 14. Each 
point is an average over M realizations, with M = 90 for 
N = 1000, M = 50 for N = 2000, M = 40 for TV = 4000, 
M = 30 for N = 8000, M = 6 for N = 16000, M = 4 
for N = 32000, and M = 4 for N = 10 5 . The error bars 
indicate the corresponding estimated error on the mean. 
The continuous line (oc TV -0 - 4 ) is a best-fit to the data, 
while the dashed lines corresponds to iV -1 / 2 and TV -1 / 3 . 


study the stability of our “final” measures of the different 
quantities, and to check for the possible role of finite N 
effects. 

Shown in Fig. [T3]is the result for the measured angular 
momentum L%, at the indicated times, in a series of sim¬ 
ulations for the different indicated values of N. As antici¬ 
pated, in the cases where different (and longer) times have 
been analyzed the measured angular momentum shows no 
significant dependence on the time: after the few dynami¬ 
cal time there is negligible interaction between the ejected 
and bound mass and the angular momentum of each are 
constant. As a function of N, on the other hand, there 
is a clear monotonic decrease, well fitted by a power law 
dependence, oc iV -7 , with 7 « 0.4 The dashed lines in¬ 
dicate, for comparison, the behaviors iV -1 / 2 and TV -1 / 3 . 

Shown in Fig. [13] is the result for the same data as 
in the previous figure, and for the different N shown, for 
the flatness ratio iso- In contrast we see in this case that 
the value measured, for all but the largest value of N, 
shows visible dependence on time. More specifically the 
bound mass evolves towards a more spherical distribution 
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in time, at a rate which is clearly faster for the smaller TV. 
Such an TV dependence of the time scale of the evolution 
of the virialized structure is clearly qualitatively coher¬ 
ent with collisional relaxation, and indeed a quantitative 
ana l ysis of such behavior from similar initial conditions in 
Theis & Spurzem ( 1999h shows that the timescale of this 


evolution is indeed in agreement with that predicted by 
two-body collisionality. At the very largest TV, on the other 
hand, the result appears to converge to a time independent 
value. In the small range of (larger) TV in which this time- 
dependence due to two body collisionality is small the fi¬ 
nal value appears to be consistent with an TV-independent 
behavior, but clearly extrapolation to larger TV would be 
needed to reach a firmer conclusion. We underline that, 
in contrast, we do not have this problem for the determi¬ 
nation above of the scaling with TV of L b p as this quantity 
is essentially frozen, as we have seen, at about t ss 2, and 
thus not modified b y the two body col li sion a lity at longer 
times. We note that I Joyce. Marcos Sz Svlos Labiiii ( 2009h 
report tests for collisionality during the collapse phase, for 
the very extreme case of a fla t (a = 0 ) profile, and con¬ 
clude, in line with the study of lAarseth. Lin fc Papaloizou 
( 1988h . that for TV > 10 3 collisional effects during the col¬ 
lapse phase are indeed negligible. This result, consistent 
with theoretical prediction for the collisional time scale, 
is borne out by the observation that the evolution of the 
macroscopic features of the collapse — potential energy, 
viral ratio, size of structure, and in particular the frac¬ 
tion of the mass ejected — are stable when TV increases, 
and likewise when the gravitational force softening is var¬ 
ied, provided it is sufficiently small to resolve the mean 
gravitational field. 


Given that, as we have emphasized, the breaking of 
spherical symmetry in the initial conditions is due to fi¬ 
nite TV fluctuations, the fact we find that the final angular 
momentum L v h decreases as TV increases, and in principle 
goes to zero as TV —> oo, is what one would naively ex¬ 
pect. A similar decrease in amplitude, with an exponent 
of about the same value, has been observed for the case 
a = 0 in Worrakitpoonoonl ( 201 fill . In our simulations of 
initial conditions from elliptical conditions we have not 
been able to detect, within the significant scatter of dif¬ 
ferent realizations at given TV, a clear trend for L £ to de¬ 
crease as TV increases. This is probably because of a weaker 
TV dependence in this case, as the spherical symmetry is 
“mostly” broken by t(0) ^ 0 which does not depend on TV. 
In contrast the lack of evidence for any TV dependence in 
the final i — once TV is sufficiently large so that collisional 
relaxation plays no role on the time scales probed — may 
be indicative of the true physical behavior: in the limit of 
exact spherical symmetry of an exactly cold system, it is 
expected that symmetry be broken due to so-called radial 
orbit instability. In future work we will study more fully 
the TV dependence of both the breaking of symmetry, and 
that of the angular momentum of the ejected particles, 
and their connection to one another. 



N 


Fig. 13. Flatness ratio i§o as a function of particle number 
TV, for the same simulations in the previous figure (a = 1) 
for the indicated values of TV. For this quantity the ob¬ 
served TV dependence is due to two body relaxation which 
make the system evolve in time towards a more spherical 
quasi-equilibrium. 


4. Discussion and conclusions 

The origin of the angular momentum in astrophysical 
structures — notably that measured observationally in 
galaxies - is a fascinating open problem in astrophysics 
and cosmology. In this article we have discussed and stud¬ 
ied a way, different to commonly considered ones such as 
tidal torques, in which angular momentum can be gener¬ 
ated by self-gravity only. It can occur in principle for initial 
conditions which evolve sufficiently violently — with large 
variations of the gravitational mean field — so that mass 
is ejected during relaxation towards a virialized state. We 
have shown using numerical simulations for two simple 
broad classes of initial conditions of this kind that signif¬ 
icant angular momentum can indeed be generated in this 
way, orders of magnitude larger in most cases that the 
angular momentum error due to numerics. 

We have, for illustrative purposes, and for simplic¬ 
ity, considered initial conditions very close to spherical 
symmetry. Indeed in all our initial conditions rotational 
symmetry is only broken fully by the Poissonian density 
fluctuations associated with the finite number of parti¬ 
cles. Interestingly we find that, despite the “artificiality” 
of our initial conditions, which are not motivated by any 
detailed specific physical model, we obtain angular mo¬ 
menta which are of comparable size — only smaller by a 
factor of two in some cases — than those typically esti¬ 
mated for galaxies. Indeed we have obtained values of the 
normalized “spin parameter” up to A ~ 0.02, while, for ex- 
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ample, Hernandez et al.l ( 2007 1 have estimated its average 
value in the 11597 spirals and elliptical galaxies observed 
by the Sloan Digital Sky Survey to be Ao = 0.04 ± 0.005, 
(and standard deviation o\ = 0.51 ± 0.05). While obvi¬ 
ously galaxy formation in particular involves much more 
complicated non-gravitational processes than described by 
the simple models we have considered here, and further 
our specific initial conditions are quite ad hoc (and not 
directly motivated by a cosmological model) it is intrigu¬ 
ing that the order of magnitude of the values we obtain 
are in line with those observed. 

We note that studies of dark m atter halos in co s molog ¬ 


ical N-body s i mulat i ons (see e.g . IVitvitska et al.l (|2002l) : 


g. V ltvit 

Bullock et al.l (2001); iBctt et all ( 2007l P give rise to sim¬ 
ilar values of spin parameter (A ~ 0.04), i.e. of the same 
order as those observed in elliptical galaxies, and as those 
we find in our simulations in certain cases. As dark mat¬ 
ter halos of standard cold dark matter cosmologies form 
through a hierarchical process rather than in the kind of 
monolithic collapse we have studied, the mechanism we 
have studied has likely no relevance in this particular con¬ 
text. Indeed it is the strong violence of the collapse leading 
to mass ejection which is crucial, and the formation of ha¬ 
los in a cold dark matter cosmolog y setting is not of thi s 
kind. Neverthe less, as discussed in ICarucci et al. ( 2014 1: 
ISamsimil ( 2015lb significant mass ejection can occur in a 
cosmological setting when halos merge and it would be 
interesting to investigate whether this could lead also to 
generation of angular momentum at a significant level, 
compared to the processes of mass accretion, f or ex am- 
ple, which has been argued in Vitvitska et al.l (2002]) to 
account for the spin of dark matter halos. 

Just as in cosmological simulations we underline that 
we are simulating essentially the collisionless regime of 
the gravitational dynamics (apart from the two body col- 
lisional effects which we observed to be present at longer 
times in smaller N simulations). In so far as this is true, 
the particle number N is relevant, in principle, in the prop¬ 
erties of the final state just because it fixes the amplitude 
of the initial fluctuations. Indeed these fluctuations are, 
for our initial conditions in this paper, simply Poissonian, 
and thus their amplitude at any scale (and all their sta¬ 
tistical properties) are regulated by this single parameter. 
Alternatively one could consider setting up an initial par¬ 
ticle with the same average density profile but with sta¬ 
tistical fluctuations with non-trivial correlation, described 
for example by non-trivial two-point correlation proper¬ 
ties. In this case one could then in principle vary N while 
keeping the relevant fluctuations fixed and obtain results 
independent of N. As the fluctuations play an essential 
role in the symmetry breaking, we expect that the detailed 
nature of such initial fluctuations may have an impact on 
the properties of the final structure. We postpone a de¬ 
tailed study of this interesting issue to future work. We 
will also explore whether cold but more irregular/clumpy 
initial conditions can produce values of the spin of the or¬ 
der of magnitude to those observed for galaxies. Further 
we will explore whether the kind of correlation we have 


discussed briefly in the last part of the paper — between 
the shape of the asymmetric virialized structure and its 
angular momentum — could be used to find observational 
evidence for or against the role of very violent relaxation 
in generating the structure of galaxies. 

Numerical simulations have been run on the Cineca 
PLX cluster (project ISCRA BSS-GC), and on the HPC 
resources of The Institute for Scientific Computing and 
Simulation financed by Region lie de France and the 
project Equip@Meso (reference ANR-10-EQPX- 29-01) 
overseen by the French National Research Agency (ANR) 
as part of the Investissements d’Avenir program. 
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